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NOMENCLATURE 


a^  distance  between  midchord  and  elastic  axis  in  semichords,  positive 

toward  the  trailing  edge 
c  full  chord  length  of  airfoil 

Cj  lift  coefficient 

Cm  moment  coefficient  about  leading  edge 

Cjtot  slope  of  the  lift  coefficient  curve 
slope  of  the  moment  coefficient  curve 
GJ  torsional  rigidity  of  wing  section 

kc  reduced  frequency  defined  as  wc/U 

k  (y+1 )M2 ,  transonic  flow  parameter 

free  stream  Mach  number 
p  ratio  of  the  semi  span  to  chord 

q  dynamic  pressure  defined  as  1/2  pU2 

l)  free  stream  velocity 

a  angle  of  attack 

0  1/O-M^2)^,  Prandtl-Glauert  Number 

Y  ratio  of  specific  heats 

6xs  shock  movement  in  chord  per  unit  perturbation 
e  perturbation  parameter 

p  free  stream  air  density 

x  ratio  of  the  maximum  airfoil  thickness  to  chord 

<}>  disturbance  velocity  potential 

(o  circular  frequency  of  oscillation 

(0)  denote  parameters  corresponding  to  base  flow 

(1)  denote  parameters  corresponding  to  calibration  flow 
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SECTION  I 
INTRODUCTION 


In  recent  years  there  is  an  increasing  trend  that  aircraft  be  operated 
at  a  speed  in  high  subsonic  or  transonic  regimes.  In  transonic  flight, 
smal 1-ampl itude  oscillations  of  a  wing  can  produce  large  variations  in  the 
aerodynamic  forces  and  moments  acting  on  that  wing.  Moreover,  phase  differences 
between  the  motion  and  the  resulting  forces  and  moments  can  be  large.  These 
characteristics  tend  to  increase  the  probability  of  encountering  aeroelastic 
instabilities.  Thus,  the  transonic  regime  has  become  a  sensitive  one  for  aero¬ 
elastic  analysis.  Studies  in  transonic  aeroelasticity  have  begun  recently. 

In  References  1  to  5,  Yang  et  al.  have  conducted  various  aeroelastic 
studies  of  airfoils  in  two  dimensional  transonic  flows.  In  Reference  6, 

Eastep  and  Olsen  have  conducted  the  flutter  analysis  of  a  rectangular  wing 
in  three  dimensional  transonic  flows.  A  state-of-the-art  review  on  the 
aeroelastic  applications  of  transonic  aerodynamics  is  given  by  Ashley  in 
Reference  7. 

In  transonic  aeroelastic  studies, the  computation  of  aerodynamic  forces  is 
a  major  task.  Previous  transonic  aeroelastic  studies  have  shown  that  the 
computer  time  required  to  compute  the  aerodynamic  data  is  quite  high.  It  is 
due  to  the  fact  that  the  transonic  aerodynamics  depends  on  flow  parameters  such 
as  Mach  number,  angle  of  attack,  airfoil  configuration,  reduced  frequency,  etc. 
in  a  nonlinear  fashion.  Because  of  this,  the  aerodynamic  computation  has  to  be 
repeated  whenever  any  one  of  the  above  flow  parameters  is  varied.  Also 
computations  become  more  complicated  when  shocks  are  present. 
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The  linear  subsonic  and  supersonic  aerodynamic  equations  yield  closed 
form  relations  between  aerodynamic  forces  and  flow  parameters.  On  the  other 
hand,  nonlinear  transonic  aerodynamic  equations  do  not  directly  yield  any 
such  relations.  Lack  of  such  closed  form  relations  has  put  restriction  on 
certain  transonic  aeroelastic  studies.  For  example,  in  order  to  compute  the 
transonic  divergence  speed  of  a  wing,  it  is  required  to  know  the  aerodynamic 
forces  as  a  function  of  angle  of  attack. 

Recently  some  techniques  based  on  perturbation  method  have  been  developed 
to  avoid  the  repetition  of  the  aerodynamic  computation  when  flow  parameters 
are  varied.  Such  techniques  have  also  led  to  simple  relations  between  aero¬ 
dynamic  forces  and  flow  parameters.  One  such  technique  is  to  use  the  concept 
of  the  strained  coordinate  perturbation  procedure  in  obtaining  the  transonic 
aerodynamic  loads. 

The  basic  concept  of  the  method  of  strained  coordinate  perturbation  is  to 
minimize  the  actual  number  of  separate  calculations  required  in  a  particular 
application  by  extending,  over  some  parametric  range,  the  usefulness  of  each 
individual  solution  determined  by  some  computational ly-expensive  procedure. 
Coordinate  straining  introduced  into  transonic  aerodynamics  as  a  means  of 
accounting  properly  for  the  movement  of  discontinuities  (shock  waves)  due  to 
changes  in  some  geometric  or  flow  parameter  is  shown  to  result  in  an  accurate 
perturbation  predictions  in  the  vicinity  of  the  discontinuity.  Detailed  studies 
of  using  this  method  as  an  effective  tool  for  reducing  the  computational 
requirements  in  transonic  aerodynamics  has  begun. 
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The  basic  concepts  of  coordinate  straining  is  given,  for  example,  in 
References  8  and  9.  An  evaluation  of  the  strained  coordinate  perturbation 
procedure  as  applied  to  nonlinear  subsonic  and  transonic  flows  is  given  by 
Stahara  et  al-in  Reference  10. 

In  Reference  10,  a  procedure  of  defining  a  unit  perturbation  by  employing 
two  nonlinear  solutions  which  differ  from  one  another  by  a  nominal  change 
in  some  geometric  or  flow  parameter,  and  then  using  that  unit  perturbation  to 
predict  a  family  of  related  nonlinear  solutions  over  a  range  of  parameter 
variation  is  discussed.  Coordinate  straining  is  used  in  determining  the  unit 
perturbation  to  account  for  movement  of  shocks  due  to  the  perturbation.  Based 
on  full  potential  solutions,  perturbation  results  are  presented  for  flows 
past  both  isolated  airfoils  and  compressor  cascades  involving  a  variety  of 
flow  and  geometry  parameter  changes  in  transonic  regime.  Comparisons  with 
corresponding  "exact"  nonlinear  solutions  indicate  a  good  accuracy  and  range 
of  validity  of  such  a  method. 

In  Reference  11,  Nixon  illustrated  a  procedure  of  perturbing  a  transonic 
flow  with  shock  wave.  The  method  is  based  on  the  use  of  a  distorted  airfoil 
as  the  initial  case  rather  than  the  real  physical  airfoil.  The  distortion  is 
chosen  such  that  the  shock  location  is  unchanged  by  the  perturbation.  The 
distorted  airfoil  is  obtained  by  the  use  of  a  strained  coordinate  system.  The 
procedure  yielded  an  algebraic  similarity  relation  between  related  airfoils 
with  shock  waves  at  different  locations.  Results  are  illustrated  for  a  NACA 
0012  airfoil  and  10%  thick  parabolic  arc  at  transonic  Mach  numbers.  The 
pressure  distributions  around  the  perturbed  airfoils  computed  by  using  both 
the  extended  integral  equation  method  and  the  perturbation  method  compare  well. 
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In  Reference  12,  Nixon  extended  the  fundamental  concept  of  Reference  11 
to  two-dimensional  lifting  flows  and  three-dimensional  lifting  flows  with 
multiple,  intersecting  shock  waves.  The  application  of  strained  coordinate 
perturbation  method  in  computing  the  steady  transonic  aerodynamic  loads  when 
Mach  number  and  angle-of-attack  are  varied  is  illustrated.  From  this  method, 
transonic  steady-state  aerodynamic  solution  at  any  given  value  of  a  parameter 
can  be  expressed  as  a  function  of  the  solutions  obtained  at  two  different 
values  of  the  same  parameter.  The  solutions  required  are  termed  as  base  and 
calibration  solutions,  respectively.  The  procedure  is  based  on  the  assumptions 
that  shocks  are  neither  created  nor  destroyed  within  the  range  of  variation 
of  the  parameter  (Mach  number  or  angle-of-attack)  and  the  perturbation  of  the 
parameter  is  small.  Using  the  method,  simple  algebraic  expressions  were  derived 
for  lift,  pitching  moment  and  drag  coefficients  in  transonic  regime.  The  method 
was  illustrated  for  a  10%  thick  parabolic  arc  and  a  NACA  64A410  airfoil  at 
transonic  Mach  numbers. 

In  Reference  13,  Nixon  compared  the  strained  coordinate  interpolation  method 
used  for  transonic  flow  in  Reference  12  with  normal  interpolation/extrapolation 
procedures.  It  was  found  that  both  methods  are  essentially  equivalent  in 
smooth  regions  of  the  solution.  However,  the  normal  linear  extrapolation  may 
not  be  applicable  in  the  region  just  behind  the  shock  wave.  The  strained 
coordinate  method  does  move  the  shock  and  its  associated  shock  foot  singularity 
to  the  correct  location  and  scales  the  strength  of  the  singularity  according  to 
linear  interpolation. 
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Based  on  the  strained  coordinate  perturbation  procedure  the  computational 
time  required  for  transonic  aeroelastic  studies  can  considerably  be  reduced 
when  parameters  such  as  Mach  number,  angle-of-attack,  thickness,  etc.  are 
varied.  Using  the  developments,  for  example,  in  Reference  12,  the  repetition 
of  steady  state  aerodynamic  computations  can  be  avoided  when  flow  parameters 
are  varied.  Thus,  there  will  be  considerable  reduction  in  the  computational 
time  for  aeroelastic  studies  using  the  steady  state  computations.  For 
example,  in  computing  static  divergence  speed  of  a  wing  only  steady  state 
aerodynamic  data  is  required. 

Transonic  unsteady  computations  based  on  time  integration,  indicial 
and  harmonic  methods  require  fairly  accurate  steady  state  initial  conditions. 
It  is  noted  that  in  both  indicial  and  harmonic  methods, the  unsteady  solution 
is  treated  as  a  small  linear  perturbation  about  a  nonlinear  steady  state 
solution.  When  flow  parameters  are  varied  in  unsteady  computations, the 
required  steady  state  initial  conditions  can  be  economically  obtained  by 
using  the  strained  coordinate  perturbation  method. 

Certain  aeroelastic  computations  require  the  aerodynamic  forces  as  a 
function  of  flow  parameters.  For  example,  in  computing  static  divergence 
speed  of  a  wing,  it  is  necessary  to  know  the  aerodynamic  forces  as  a  function 
of  angle-of-attack.  Such  functions  for  transonic  regime  can  be  obtained  by 
the  strained  coordinate  perturbation  method. 

In  this  study,  some  preliminary  applications  of  the  strained  coordinate 
perturbation  procedure  in  transonic  aeroelasticity  are  investigated.  The 
main  emphasis  is  on  reducing  the  computational  time  in  transonic  aeroelastic 
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studies.  The  developments  In  Reference  12  are  used  for  this  purpose. 

First,  the  use  of  the  strained  coordinate  perturbation  method  in  computing 
the  steady  and  unsteady  aerodynamic  pressure  distributions  for  varying  flow 
parameters  were  studied.  As  a  computational  example,  a  parabolic  arc  was 
selected.  Variations  of  angle-of -attack,  Mach  number  and  thickness  were 
considered.  The  aerodynamic  solutions  at  base  and  calibration  values  of  the 
flow  parameter  were  computed  by  using  the  transonic  code  LTRAN2  developed  by 
Ballhaus  and  Goorjian  (Reference  14). 

For  a  6%  parabolic  arc  at  Mach  number  0.8,  the  steady  and  unsteady 
aerodynamic  pressure  coefficients  were  computed  at  angle-of-attack  0.8°  by 
strained  coordinate  perturbation  method.  This  was  based  on  the  base  and 
calibration  angles-of-attack  equal  to  0.4°  and  0.6°,  respectively.  The  aero¬ 
dynamic  results  obtained  at  angle-of-attack  0.8°  by  strained  coordinate 
perturbation  method  were  compared  with  those  directly  obtained  by  LTRAN2. 

The  comparison  is  good. 

For  the  6%  parabolic  arc  at  zero  angle-of-attack,  the  steady  and  unsteady 
aerodynamic  results  were  computed  at  Mach  number  0.86.  This  was  based  on 
the  base  and  calibration  Mach  numbers  equal  to  0.854  and  0.856,  respectively. 
Present  results  are  compared  with  those  directly  obtained  by  LTRAN2.  The 
comparison  is  good. 

For  a  6.5%  thick  parabolic  arc  at  Mach  number  0.85,  aerodynamic  pressure 
coefficients  were  computed  by  strained  coordinate  method.  The  base  and 
calibration  thicknesses  were  6%  and  7%,  respectively.  Present  results  compare 
well  with  those  obtained  directly  by  LTRAN2. 
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Finally,  the  use  of  the  strained  coordinate  perturbation  method  in 
computing  the  transonic  divergence  dynamic  pressure  of  a  typical  slender 
straight  wing  with  conventional  airfoil  is  illustrated.  Since  two-dimensional 
aerodynamics  were  used  in  obtaining  the  aerodynamic  loads,  the  aeroelastic 
equations  were  derived  by  using  strip  theory.  The  solution  for  divergence 
dynamic  pressure  was  obtained  by  Rayleigh's  method  (Reference  15). 

As  a  computational  example,  a  10%  thick  parabolic  arc  airfoil  section 
was  selected.  By  using  the  aerodynamic  force  coefficients  computed  by  the 
strained  coordinate  perturbation  method,  the  effect  of  Mach  number  on  transonic 
divergence  dynamic  pressure  was  studied.  These  results  were  compared  with 
those  obtained  by  linear  subsonic  aerodynamic  theory. 

Based  on  the  present  procedure  a  computer  program  to  post  process 
LTRAN2  aerodynamic  data  for  strained  coordinate  perturbation  method  was 
written.  The  listing  with  user's  manual  is  given.  This  program  creates  new 
steady  state  initial  conditions  for  LTRAN2  by  using  the  base  and  calibration 
steady  state  solutions  obtained  by  LTRAN2  for  varying  the  flow  parameters. 
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SECTION  II 


STRAINED  COORDINATE  PERTURBATION  EQUATIONS 


In  this  section, the  strained  coordinate  perturbation  equations  are 
presented.  The  equations  are  based  on  the  assumptions  that  shocks  are 
neither  created  nor  destroyed  during  perturbation  and  the  order  ,f  the 
perturbation  is  small.  The  main  objective  is  to  obtain  from  two  or  more 
solutions  an  algebraic  relation  that  connects  the  flow  variables  for  a  range 
of  one  or  more  parameters,  thus  leading  to  a  rapid  computation  of  these 
related  flows.  The  effect  of  shock  movement  during  perturbation  is  accounted 
by  using  the  procedure  given  in  Reference  12. 

Equations  are  presented  for  two-dimensional  transonic  inviscid  flow 
governed  by  small -disturbance  conditions. 

The  basic  steady-state  equation  in  a  scaled  form  is  given  by 


where  (x,z)  is  a  Cartesian  Coordinate  System,  with  x  aligned  with  the  airfoil 
chord  and  related  to  the  physical  coordinate  system  (x,z)  by  the  transformation 

x  =  x,  z  =  6 z  (2) 

where,  if  is  the  freestream  Mach  number,  then 

6  =  1.0/(1  -  M*)1/2  (3) 

The  potential  <j>  (x,z)  is  expanded  as  a  series  in  small  parameter  e  such 
as 

<p  (x,z)  =  <P0  (x',z)  +  e  <Pl  (x',z)  + -  (4) 
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where  x'  is  the  strained  coordinate.  Shock  is  assuifed  to  be  normal  to  the 
freestream,  thus  only  x  coordinate  straining  is  required.  Following  the 
discussion  given  in  Reference  12  the  strained  coordinate  system  is  defined  by 


x  =  x‘  +  t  6x$xt  (x1 ) 


(5) 


where,  if  x^  is  the  location  of  the  shock  in  the  (x',z)  coordinates,  then 


x.  ( x ’ )  =  x1  (1  -  x1 

x  , . 
s  s 


( 1  -  x ’ )  ;  0<  x'  < 

H  -  x'j 


(6a) 


x^  ( x ’ )  =  0;  x’>l,  x'<0  (6b) 

and  efix$  is  the  amount  by  which  the  shock  moves  during  the  perturbation. 

If  by  some  method,  say  finite  difference,  two  solutions  termed  as  base 
and  calibration  are  known,  then  an  expression  for  the  disturbance  velocity 

potential  can  be  expressed  as 

(o) 

♦x(x,z)  =  <J>X ,  (x'.z)  [1  -  e6xsxlx,  (x1 )]  + 


e4>lx.  (x1  ,z ) 


(7) 


where 


(1)  „  (o) 

e0*ix(x\z)  =  <l>x  (x,z)  -  <>x,  (x'.z)  [1  -  eo6xsxlx,(x')] 


X 

X 


c  fix 

0 


s 


♦ 


(o) 

x' 


=  X'  +  eo5xsx1(x') 

*  x'  +  efix$x^(x' ) 

*  the  change  in  the  shock  location  between  base  and 
calibration  solutions 

*  base  solution 

*  calibration  solution 
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The  expressions  for  the  perturbations  e  and  eQ  depend  on  the  flow  or 
geometry  parameter  perturbed.  Following  the  discussion  given  in  Reference  12 
for  Mach  number  variation  the  expressions  for  e  and  e0  can  be  written  as 


~(y+1)m!  .  (y+1)m2 

/ 

' (y+1 )h3 

(1-M£)3/2  (1-M*)3/2 

*1 

[o^/2\ 

(8a) 


= r (y+i  )mi 


(y+l)M 

(1-Mq) 


(y+1 )m3 
(1^)3/2 


(8b) 


where  y  =  ratio  of  specific  heats 

M0  =  freestream  Mach  number  of  base  flow 

Mj  =  freestream  Mach  number  of  calibration  flow 

M^  =  current  freestream  Mach  number 

q  =  the  exponent  in  the  transonic  parameter  taken  to  be  2 
For  angle-of -attack  variation  the  expressions  for  e  and  c0  can  be  written 


as 


£  =  a  -  aQ 


-  a. 


(9a) 

(9b) 


where 

a  =  current  angle-of-attack 
cij  =  calibration  angle-of-attack 
a  =  base  angle-of-attack. 
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For  thickness  variation  the  expressions  c  and  c0  can  be  written  as 
c  =  t  -  t0  (10a) 

e0  s  ti  -  x0  (10b) 

where  t  »  current  ratio  of  maximum  thickness  to  chord  of  airfoil 


t i  s  ratio  of  maximum  thickness  to  chord  for  calibration  airfoil 
t0  *  ratio  of  maximum  thickness  to  chord  for  base  airfoil. 

Based  on  the  Equation  7  expressions  for  lift  and  moment  coefficients 
can  be  written  as 


(11a) 


(lib) 


where  C*,  Clo  and  Cil  are  lift  coefficients  corresponding  to  current,  base 
and  calibration  values  of  the  flow  parameters,  respectively;  C„,,  CmQ  and 
C,^  are  the  moment  coefficients  corresponding  to  current,  base  and 
calibration  values  of  the  flow  parameters,  respectively;  the  transonic  flow 
parameter  k  is  defined  as  (y+1)  fl*,  and  (o)  and  (1)  denote  computations 
made  at  base  and  calibration  values  of  the  flow  parameters, respectively. 


11 


SECTION  III 


APPLICATION  TO  STEADY  AND  UNSTEADY  COMPUTATIONS 

It  is  a  common  practice  to  conduct  parametric  studies  in  transonic 
aeroelasticity.  Usually  flow  parameters  such  as  Mach  number  and  angle 
of  attack  are  varied.  In  such  cases  the  repetition  of  steady  state 
computations  can  be  avoided  by  using  the  strained  coordinate  perturbation 
method.  Thus,  there  will  be  direct  saving  in  the  computational  time  for 
aeroelastic  studies  that  depend  on  the  steady  state  computations. 

The  unsteady  computations  based  on  time  integration,  indicial  and 
harmonic  methods  depend  on  accurate  steady  state  initial  conditions.  For 
given  values  of  flow  parameters,  the  steady  state  computation  requires  about 
25  to  50  percent  of  computer  time  in  transonic  computer  codes  such  as 
LTRAN2.  Considerable  amount  of  computer  time  can  be  saved  by  computing  the 
steady  state  initial  conditions  from  the  strained  coordinate  perturbation 
method. 

In  this  section,  results  are  shown  to  illustrate  the  use  of  the  strained 
coordinate  perturbation  method  in  steady  and  unsteady  computations  when  Mach 
number,  angle-of-attack  and  thickness  ratio  are  varied.  The  airfoil  considered 
is  a  parabolic  arc.  The  Mach  numbers  considered  were  in  the  transonic 
range.  The  base  and  calibration  solutions  were  obtained  from  LTRAN2  by 
using  a  finite  difference  mesh  with  119  horizontal  grid  points  and  79 
vertical  grid  points, 
a.  Mach  Number  Variation 

A  parabolic  airfoil  of  6%  maximum  thickness-chord  ratio  at  zero  angle 
of  attack  was  considered.  Based  on  the  base  and  calibration  Mach  numbers 
equal  to  0.854  and  0.856,  respectively,  results  for  Mach  number  0.86  was 
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obtained  by  strained  coordinate  perturbation  method.  They  are  compared 
with  those  obtained  directly  from  LTRAN2. 

Figure  1  shows  the  steady  state  pressure  coefficient  curve  obtained 
by  strained  coordiante  perturbation  method  for  M  =  0.86  by  using  Equation 
7.  The  values  for  e  and  e0  were  obtained  by  Equations  8a  and  8b  as  0.0748 
and  0.0240,  respectively.  The  amount  of  shock  movement  e06xs  between  the 
base  and  calibration  Mach  numbers  was  equal  to  0.02  chord.  The  computational 
time  required  to  obtain  this  curve  was  about  5  seconds  on  CYBER  74  computer. 

In  Figure  1,  the  steady  state  pressure  curve  obtained  directly  from  LTRAN2 
is  also  shown.  The  curve  obtained  by  strained  coordinate  perturbation  method 
agreeswell  with  that  obtained  by  LTRAN2.  There  are  some  discrepancies  near 
the  shock.  They  may  be  due  to  the  lack  of  fine  enough  grid  near  the  shock. 

Based  on  the  base  and  calibration  steady  state  solutions  obtained  at 
Mach  numbers  0.854  and  0.856,  respectively,  initial  conditions  for  Mach 
number  0.86  were  computed  by  the  strained  coordinate  perturbation  Equation  7. 
This  was  carried-out  by  using  a  computer  program  given  in  the  Appendix.  Using 
this  steady  state  initial  condition,  unsteady  results  were  computed  by  the 
time  integration  procedure  of  LTRAN2.  The  reduced  frequency  kc  was  assumed 
as  0.1. 

Figure  2  shows  two  sets  of  unsteady  pressure  curves  obtained  by  using  the 
present  initial  condition  and  the  initial  condition  directly  obtained  by  LTRAN2. 
The  curves  were  obtained  at  non-dimensional  time  ut  equal  to  18.06  radians, 
where  ui  is  the  circular  frequency  of  oscillation. 

In  Figure  2,  it  can  be  observed  that  results  obtained  by  the  two  methods 
agree  well.  Slight  discrepancies  in  the  initial  steady  state  conditions 
obtained  by  the  present  method  do  not  have  any  influence  on  the  unsteady  pressure 
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Figure  1.  Effect  of  Mach  Number  Perturbation  for  6%  Thick  Parabolic 
Arc  on  Steady  Pressure  Curves  at  M  =  0.86. 


UNSTEADY  PRESSURE  COEFFICIENT 
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Figure  2.  Effect  of  Mach  Number  Perturbation  for  6%  Thick  Parabolic 
Arc  on  Unsteady  Pressure  Curves  at  M  =  0.86. 


curve.  Similarly,  unsteady  lift  and  moment  coefficients  compared  well 
between  the  two  methods. 

It  was  noticed  that  to  obtain  unsteady  results  for  3  cycles,  the 
present  method  required  about  60%  of  the  computer  time  required  for  direct 
LTRAN2  computations.  Based  on  the  same  base  and  calibration  Mach  numbers, 
steady  and  unsteady  results  can  be  obtained  for  Mach  numbers,  say,  between 
0.85  and  0.86.  Thus,  there  can  be  a  considerable  saving  in  the  computer  time 
if  additional  aeroelastic  computations  have  to  be  conducted  in  this  range, 
b.  Angle  -of-Attack  Variation 

The  6%  thick  parabolic  arc  at  M  =  0.8  was  considered.  Steady  and 
unsteady  aerodynamic  results  for  angle-of-attack  0.8°  was  obtained  by 
strained  coordinate  perturbation  procedure.  The  base  and  calibration  angles- 
of-attack  were  equal  to  0.4°  and  0.6°,  respecti vely.  Present  results  are 
compared  with  those  directly  obtained  by  LTRAN2. 

The  values  for  e  and  e0  were  obtained  from  Equations  9a  and  9b  as  0.4 
and  0.2,  respectively.  Since  there  is  no  shock,  the  value  for  e06xs  is 
equal  to  zero. 

Figure  3  shows  the  steady  pressure  curves  obtained  by  both  present 
procedure  and  LTRAN2.  The  agreement  between  the  methods  is  excellent.  The 
computational  time  required  for  the  present  procedure  was  about  5  seconds  on  the 
CYBER  74  computer. 

Based  on  the  base  and  calibration  solutions,  steady  state  initial 
conditions  for  unsteady  computation  was  obtained  by  Equation  7  for  angle- 
of-attack  0.8°.  This  was  carried-out  by  using  the  computer  program  given  in  the 
Appendix.  Using  this  initial  condition,  unsteady  computations  were  conducted 
at  reduced  frequency  kc  =  0.1  by  pitching  the  airfoil  at  the  midchord. 
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STEADY  PRESSURE  CURVE 
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X/C 

Figure  3.  Effect  of  Angle  of  Attack  Perturbation  for  6%  Thick 
Parabolic  Arc  on  Steady  Pressure  Curves  at  a  =  0.8°. 
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From  unsteady  computations  it  was  found  that  the  present  method  gives 
results  identical  to  those  obtained  by  LTRAN2  directly.  This  is  because 
of  the  excellent  comparison  obtained  between  the  two  corresponding  steady 
state  curves. 

It  was  noticed  that  to  obtain  unsteady  results  for  3  cycles,  the  present 
method  required  about  40%  of  the  computer  time  required  for  direct  LTRAN2 
computations.  Based  on  the  same  base  and  calibration  angles-of-attack,  steady 
and  unsteady  results  can  be  obtained  for  angles-of-attack,  say,  between  0.0° 
and  1.0°  from  the  present  procedure.  Thus,  there  can  be  considerable  saving 
in  the  computer  time  if  aeroelastic  computations  have  to  be  conducted  in 
this  range. 

c.  Thickness  Variation 

A  parabolic  arc  airfoil  at  Mach  number  0.85  with  zero  angle-of-attack 
was  considered.  Based  on  base  and  calibration  maximum  thickness  to  chord 
ratios  equal  to  0.06  and  0.07,  respectively,  results  for  maximum  thickness 
to  chord  ratio  equal  to  0.065  was  obtained  by  strained  coordinate  perturbation 
method.  They  are  compared  with  those  obtained  directly  from  LTRAN2. 

Figure  4  shows  the  steady  state  pressure  coefficient  curve  obtained  by 
strained  coordinate  perturbation  method  for  x  =  0.065  by  using  Equation  7. 

The  values  for  e  and  e0  were  obtained  from  Equations  10a  and  10b  as  0.005 
and  0.01,  respectively.  The  amount  of  shock  movement  e0<5xs  between  the  base 
and  calibration  thickness  ratios  was  equal  to  0.08  chord.  The  computational 
time  required  to  obtain  this  curve  was  about  5  seconds  on  CYBER  74  computer. 
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In  Figure  4, steady  state  pressure  curve  obtained  directly  from  LTRAN2 
is  also  shown.  The  agreement  between  two  curves  is  good.  There  is  some 
discrepancy  near  the  shock.  It  may  be  due  to  lack  of  a  fine  grid  near  the 
shock. 

Based  on  the  base  and  calibration  steady  state  solutions  obtained  at 
thickness  ratios  0.06  and  0.07,  respectively,  initial  conditions  for  thickness 
ratio  equal  to  0.065  was  computed  by  the  strained  coordinate  perturbation 
Equation  7.  This  was  carried  out  by  using  a  computer  program  given  in  the 
Appendix.  Using  this  steady  state  initial  condition,  unsteady  results  were 
computed  by  time  integration  procedure  of  LTRAN2.  The  reduced  frequency  kc 
was  assumed  as  0.1. 

Figure  5  shows  two  sets  of  unsteady  pressure  curves  obtained  by  using 
present  initial  conditions  and  initial  conditions  directly  obtained  by  LTRAN2. 
The  curves  were  plotted  at  non-dimensional  time  iot  equal  to  18.06  radians. 

In  Figure  5,  it  can  be  observed  that  results  obtained  by  the  two  methods 
agree  well.  Slight  discrepancies  in  the  initial  steady  state  conditions 
obtained  by  present  method  do  not  have  any  influence  on  unsteady  pressure  curve. 
Similarly,  the  unsteady  lift  and  moment  coefficients  compared  well  between  the 
two  methods. 

It  was  noticed  that  to  obtain  unsteady  results  for  3  cycles,  the  present 
method  required  about  60%  of  the  computer  time  required  for  the  direct  LTRAN2 
computations.  Based  on  the  same  base  and  calibration  thickness  ratios,  steady 
and  unsteady  results  can  be  obtained  for  thickness  ratios,  say,  between  0.05 
and  0.08  from  present  procedure.  Thus,  there  can  be  considerable  saving  in 
the  computer  time  if  additional  aeroelastic  computations  have  to  be  conducted 
in  this  range. 
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UNSTEADY  PRESSURE  COEFFICIENT 


0.065.  M  =  0.85.  «  =  0.0°.  Ke  =o.l 
-  PITCHING  ABOUT  MIDCHORD 


Effect  of  Thickness  Perturbation  for  6.5%  Thick  Parabolic 
Arc  on  Unsteady  Pressure  Curves  at  M  =  0.85  and  a  =  0.0°. 


SECTION  IV 


APPLICATION  TO  THE 
TRANSONIC  DIVERGENCE  PROBLEM 

Based  on  the  strained  coordinate  perturbation  method,  computational 
time  required  for  transonic  aeroelastic  studies  can  considerably  be  reduced. 

One  such  study  where  this  method  can  be  used  is  in  computing  the  transonic 
static  divergence  speeds  of  slender  straight  wings. 

In  this  study,  the  transonic  divergence  dynamic  pressure  of  a  typical 
slender  straight  wing  with  conventional  airfoil  is  computed  by  using  the 
method  of  strained  coordinates.  The  required  base  and  calibration  aero¬ 
dynamic  solutions  were  obtained  by  LTRAN2.  Since  two-dimensional  aero¬ 
dynamics  was  used  in  obtaining  the  aerodynamic  loads,  the  aeroelastic 
equations  were  derived  by  using  strip  theory.  The  solution  for  divergence 
dynamic  pressure  was  obtained  by  Rayleigh's  method  (Reference  15). 

As  a  computational  example,  a  10%  thick  parabolic  arc  airfoil  section 
was  selected.  First,  the  use  of  the  method  of  strained  coordinates was 
illustrated  at  a  transonic  Mach  number  of  0.8  for  varying  angle-of -attack. 

Then  the  effect  of  Mach  number  on  transonic  divergence  dynamic  pressure  was 
studied  for  various  positions  of  the  elastic  axis.  These  results  were  compared 
with  those  obtained  by  linear  subsonic  theory, 
a.  Formulation  of  Divergence  Equations 

It  is  assumed  that  the  wing  is  slender  and  straight  so  that  the  three 
dimensional  effects  of  the  aerodynamics  can  be  neglected  and  that  strip  theory 
can  be  used  in  deriving  an  expression  for  the  static  divergence  dynamic 
pressure.  Assuming  that  the  wing  torsional  deformation  pattern  is  invariable. 
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with  respect  to  the  load  distribution  and  using  Rayleigh's  method  an 
expression  for  the  static  divergence  dynamic  pressure  can  be  written  as 

“/  {%)*  d« 

%  '  ^ -  02) 

C'f  lCia  <>  *  V«  *  Cn»)f!d5 

where  qD  =  1/2  pU2,  is  the  divergence  dynamic  pressure;  p,  density  of  the  air; 
U,  flight  speed;  GJ,  torsional  rigidity  of  the  wing  section  assumed  to  be 
constant  along  span;  p,  ratio  of  the  span  to  chord;  f,  divergence  mode  shape 
expressed  as  a  function  of  £;  £,  is  the  ratio  of  the  distance  measured  in 
full  chords  from  the  root  to  a  span  station;  c,  full  chord  length  of  the 
wing  assumed  to  be  constant  along  the  span;  ,  slope  of  the  aerodynamic 
lifting  force  with  respect  to  the  angle-of-attack;  C  ,  slope  of  the  aero- 
dynamic  pitching  moment  (measured  about  the  leading  edge)  with  respect  to 
the  angle-of-attack;  and  a^,  position  of  the  elastic  axis  measured  in  semi- 
chords  from  the  midchord  (positive  towards  trailing  edge). 

For  the  subsonic  and  supersonic  flows  the  aerodynamic  equations  are 
linear  and  the  aerodynamic  force  coefficients,  C,  and  C  ,  depend  on  the 
angle-of-attack  linearly.  Hence,  Equation  12  can  be  integrated  once  an 
expression  for  f(£)  is  assumed.  On  the  other  hand,  for  the  transonic  flows, 
the  aerodynamic  equations  are  non-linear  and  the  aerodynamic  force  coefficients, 
and  Sn’  dePend  on  angle-of-attack  in  a  non-linear  fashion.  Hence,  it  is 
required  to  express  and  Cm  as  a  function  of  angle-of-attack  in  order  to 
integrate  Equation  12.  Such  expression  can  be  derived  for  transonic  flows 
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by  using  the  method  of  strained  coordinates. 

Assuming  that  perturbation  in  the  angle-of-attack  is  small  and  the  shock 
is  neither  created  nor  destroyed  within  the  range  of  the  perturbation, 
and  Cm  can  be  expressed  as  (See  Equations  11a  and  lib) 


ct  .  Ct  *  <e/e  )  (Ct  -  ct  ) 

0  1  0 


C,  =  C  +  (e/e  )  (C  -  C  ) 
m  m  o  m,  m 

o  i  o 


(13) 

(14) 


where  and  C£q  are  calibration  and  base  lift  coefficients,  respectively; 
Cm)  and  Cm0  are  calibration  and  base  moment  coefficients  (measured  about  the 
leading  edge),  respectively;  e  —  a  —  a0;  e0  =  oij  —  a0;  is  the  calibration 
angle-of-attack;  and  a0  is  the  base  angle-of-attack. 

Equations  13  and  14  lead  to  simple  expression  for  and  Cma  as 


C*a  =  (C£i  "  Ct0)/(a,  -  a0)  (15) 

=  (cmi  “  Cni0)/(al  ao)  (16) 

Since  and  Cma  given  by  Equations  15  and  16  are  independent  of  a. 
Equation  12  can  be  integrated  for  a  known  function  of  f(S)>  This  function 
should  approximately  represent  the  divergent  mode  that  satisfies  the  boundary 
conditions,  namely,  f (0)  =  0  and  f ’ (p )  *  0.  One  such  function  is  given  by 

f(0  *  2S/P  -  (S/P)2  (17) 
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Substituting  Equations  15,  16  and  17,  into  12  and  integrating  yields 

qD  *  15/6A  (18) 

where  qp  =  qp  p2  cVGJ,  non-dimensional  static  divergence  dynamic  pressure 
and  A  =  (C^  -  C^)  (1  +  ah)/2(a,  -  aQ)  +  (Cm^  -  -  a.)  (19) 

b.  Results 

The  Mach  numbers  considered  in  this  study  are  in  transonic  range  for  the  10" 
thick  parabolic  arc. 

First,  a  case  with  angle-of-attack  varying  from  0°  to  1.0°  at  Mach  number 
0.80  was  considered  to  illustrate  the  use  of  the  method  of  strained  coordinates. 
Lift  and  moment  coefficients  were  computed  by  both  the  finite  difference  and 
the  strained  coordinate  methods  and  they  are  compared. 

Figure  6  shows  the  plots  of  lift  coefficient  versus  angle-of-attack 
obtained  by  both  finite  difference  and  strained  coordinate  methods.  The  curve 
for  strained  coordinate  method  was  based  on  base  and  calibration  angles-of -attack 
equal  to  0.2°  and  0.8°,  respectively.  The  corresponding  results  for  the 
pitching  moment  coefficient  (about  the  leading  edge)  are  shown  in  Figure  7. 

In  both  Figures  6  and  7,  the  aerodynamic  computations  at  base  and  calibration 
angles-of-attack  were  also  made  by  the  successive  line  over  relaxation  method. 

Results  in  Figures  6  and  7  show  that  the  method  of  strained  coordinates 
agree  well  with  the  finite  difference  method.  The  level  of  agreement  is  better 
for  lift  coefficients  when  compared  to  that  for  moment  coefficients.  The 
total  computational  time  required  by  the  method  of  strained  coordinates  was 
about  1/5  of  that  required  for  the  finite  difference  method.  Also  the  method 
of  strained  coordinates  assumes  simple  relations  between  the  aerodynamic  force 


L  PARABOLIC  ARC  (10%)  ,  M=0-8 


Figure  6.  Variation  of  Lift  Coefficient  With  Angle  of  Attack  for  10% 
Thick  Parabolic  Arc  at  M  *  0.80. 


Figure  7.  Variation  of  Moment  (About  Leading  Edge)  Coefficient  With 
Angle  of  Attack  for  10%  Thick  Parabolic  Arc  at  M  =  0.80. 
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coefficients  and  the  angle-of-attack. 

The  application  of  the  method  of  strained  coordinates  in  computing  the 
transonic  divergence  characteristics  was  considered  next.  The  characteristics 
of  the  static  divergence  dynamic  pressures  for  varying  Mach  number  at  various 
values  of  the  position  of  elastic  axis  were  obtained.  The  base  and  the 
calibration  angles  were  assumed  as  0.2°  and  0.6°,  respectively.  The  Mach 
numbers  considered  were  0.76,  0.78,  0.79,  0.80,  0.805,  and  0.81,  respectively. 

Table  1  shows  the  lifting  force  and  pitching  moment  (about  the  leading 
edge)  coefficients  obtained  at  base  and  calibration  angles-of-attack  for  six 
Mach  numbers.  It  is  observed  in  the  table  that  both  the  coefficients  increase 
non-linearly  with  increase  in  Mach  number.  In  the  same  table  corresponding 
coefficients  obtained  by  linear  aerodynamic  theory  are  also  given  for 
reference.  The  differences  in  the  values  are  mainly  due  to  the  presence 
of  shocks. 

Based  on  the  aerodynamic  coefficients  given  in  Table  1  and  using  Equation 
18,  static  divergence  dynamic  pressures  were  computed.  The  values  for  the 
position  of  elastic  axis  were  assumed  as  0.0,  -0.1,  and  -0.2.  Figures  shows 
the  plots  of  divergence  dynamic  pressure  parameter  qo  versus  Mach  number.  In 
the  same  figure  the  corresponding  results  obtained  by  the  subsonic  theory 
are  also  shown. 

In  Figure  8,  it  is  observed  that  the  static  divergence  dynamic  pressure 
obtained  by  transonic  aerodynamics  Increased  with  the  increase  of  Mach  number. 
The  increase  is  more  rapid  at  higher  Mach  numbers.  Also  the  static  divergence 
dynamic  pressure  increases  as  the  elastic  axis  move  towards  leading  edge  from 
the  midchord.  The  curves  shift  to  the  left  as  the  elastic  axis  moves  toward 
the  leading  edge  from  the  midchord. 
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TABLE  1  AERODYNAMIC  COEFFICIENTS  FOR  10%  THICK 
PARABOLIC  ARC  AT  VARIOUS  MACH  NUMBERS 


MACH 

NUMBER 

j  Lift  Coef.  Cjj, 

* - -  - -  -  - 

Moment  Coef.  Cm 

CASE 

a  =  0.2° 

a  =  0.6° 

a  =  0.2° 

a  =  0.6° 

0.760 

1 

0.03893 

0.11710 

-.01139 

-.03426 

2 

0.01074 

0.03222 

-.00268 

-.00806 

1 

0.04308 

0.13033 

-.01314 

-.03994 

2 

0.01116 

0.03347 

-.00280 

-.00837 

0.790 

1 

0.04699 

0.14380 

-.01503 

-.04681 

2 

0.01139 

0.03416 

-.00285 

-.00854 

0.800 

1 

!  0.05477  i 

0.17158  j 

-.01945 

-.06275 

2 

0.01164  j  0.03491 

-.00291 

-.00873 

0.805 

'  ! 

I 

0.06222  1 

0.19965 

-.02382 

-.07989 

:  2  1 

j-. 

0.01177  ! 

0.03530 

-.00294 

-.00882 

l 

0.07322 

0.25171 

-.03050 

-.11251 

2 

0.01190 

0.0357 

-.00298 

-.00893 

1.  Transonic  Method 

2.  Subsonic  Method 
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DIVERGENCE  DYNAMIC  PRESSURE  PARAMETER 


Figure  8. 


Effect  of  Mach  Number  on  Static  Divergence  Dynamic 
Pressure  for  a  Straight  Slender  Wing  at  Various  Positions 
of  Elastic  Axis. 
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Figure  8  shows  that  the  results  obtained  by  subsonic  theory  do  not 
agree  with  those  obtained  by  the  present  transonic  theory,  especially  at 
higher  Mach  numbers.  Also  agreement  becomes  bad  as  the  elastic  axis  moves 
towards  the  leading  edge.  The  discrepancies  are  mainly  due  to  the  presence 
of  shocks  which  are  incorporated  only  in  the  transonic  theory.  It  is  noted 
that  the  rapid  changes  in  the  divergence  dynamic  pressures  at  the  higher 
Mach  numbers  are  due  to  the  movement  of  the  shock  towards  the  trailing  edge. 

The  increase  in  the  divergence  dynamic  pressures  with  the  increase  of 
Mach  number  can  further  be  explained  as  follows.  The  center  of  pressure 
(CP)  moves  from  the  quarter  chord  towards  midchord  with  increasing  Mach 
number.  Thus  CP  moves  towards  the  elastic  axis  which  is  located  near  mid¬ 
chord  for  this  study.  As  a  result,  the  divergence  dynamic  pressure  increases 
(See  Equation  8-40  of  Reference  15). 
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SECTION  V 


CONCLUDING  REMARKS 

Based  on  the  present  study  the  following  concluding  remarks  may  be 

made. 

(1)  The  steady  state  pressure  curves  obtained  by  using  strained  coordinate 
perturbation  method  compare  well  with  those  directly  obtained  by  LTRAN2. 

Some  small  discrepancies  obtained  near  the  shock  are  due  to  the  lack  of 

a  fine  enough  grid. 

(2)  The  computational  time  required  to  compute  the  steady  state  pressure 
curve  by  the  strained  coordinate  perturbation  method  for  a  known  base  and 
calibration  solution  is  about  5  seconds  on  CYBER  74  computer. 

(3)  Unsteady  results  based  on  the  initial  conditions  obtained  by  the  present 
method  compare  well  with  those  obtained  directly  by  LTRAN2. 

(4)  The  present  procedure  shows  about  40  to  50%  saving  in  the  computer  time 
for  typical  unsteady  computations  required  in  aeroelastic  analysis. 

(5)  Based  on  aerodynamic  forces  computed  by  the  strained  coordinate  per¬ 
turbation  equations,  transonic  divergence  speeds  of  slender  straight  wings 
can  be  computed. 

(6)  Based  on  the  present  computations  it  is  found  that  the  divergence  speed 
of  a  slender  straight  wing  with  10%  thick  parabolic  arc  section  increase  with 
increase  in  Mach  number.  On  the  other  hand,  linear  theory  predicts  a 
different  behavior. 

(7)  The  present  procedure  can  be  extended  to  three  dimensional  steady  and 
unsteady  transonic  computations  by  using  the  corresponding  developments  in 
strained  coordinate  methods. 
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APPENDIX 


A  computer  program  to  create  the  steady  state  initial  conditions  by 
using  the  strained  coordinate  perturbation  method  is  presented.  This  program 
is  compatible  for  LTRAN2.  This  program  can  create  a  new  steady  state  initial 

condition  based  on  base  and  calibration  steady  state  initial  conditions  obtained 
by  LTRAN2. 

a.  Description  of  the  INPUT 

(1)  Initial  DATA  Card  (415,  4F10.4)  one  card  to  define  following  parameters. 
Columns  Description  Variable. 


1-5 

Number  of  mesh  points  in 
vertical  direction 

L  MAX 

6-10 

Number  of  mesh  points  in 
horizontal  direction 

JMAX 

11-15 

Mesh  point  corresponding  to 
leading  edge 

JLE 

16-20 

Mesh  point  corresponding  to 
trailing  edge 

JTE 

21-30 

Distance  of  the  shock  from 
leading  edge  measured  in  chords 
for  base  flow 

XSO 

31-40 

Distance  of  the  shock  from 
leading  edge  measured  in  chords 
for  calibration  flow 

XS1 

41-50 

Value  of  the  perturbation 
parameter  e  (See  Section  11) 

EP 

51-60 

Value  of  the  perturbation 
parameter  eQ  (See  Section  II) 

EP0 

(2)  Mesh  Card  (8E10.4) 

Cards  to  define  JMAX  values  of  the  X  mesh  points. 


35 


(b)  Description  of  Logical  Files  Used 


TAPE  1 
TAPE  2 
TAPE  3 


Contains  steady  state  initial  conditions  from  LTRAN2  for 
base  flow  on  INPUT 

Contains  steady  state  initial  conditions  from  LTRAN2  for 
calibration  flow  on  INPUT 

Contains  the  steady  state  initial  conditions  to  LTRAN2 
for  the  current  flow  on  OUTPUT 
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nnnnnnnnnnnnnnnnnnnnnnnnn 


PROGRAM  MAIN (INPUT. OUTPUT » TAPE5=INPUT. TAPE6=0UTPUT. TAPE1 . 

1  TAPE2»  TAPE3*  TAPE4 ) 

DIMENSION  PU< 119).PL( 119).P(79» 119).DUMY1(119). DUMY2C119). 

1  X( 119) * TX( 119) 

PROGRAM  TO  CREATE  STEADY  STATE  IMTIAL  CONDITIONS  FOR  LTRAN2 
BY  STRAINED  COORDINATE  PERTURBATION  METHOD 
JMAX=NUMBER  OF  HORIZONTAL  GRID  POINTS  USED  LTRAN2 
LMAX=NUMBER  OF  UERTICAL  GRID  POINTS  USED  IN  LTRAN2 
JLE=HOR I ZONT AL  GRID  POINT  OF  LEADING  EDGE 
JTE=HOR I ZONT AL  GRID  POINT  OF  TRAILING  EDGE 
XSO=SHOCK  POSITION  IN  BASE  SOLUTION 
XSl=SHOCK  POSITION  IN  CALIBRATION  SOLUTION 
EP=UALUE  OF  PERTURBATION  PARAMETER  FOR  CURRENT  FLOW 
EPO=UALUE  OF  PERTURBATION  PARAMETER  FOR  CALIBRATION  FLOW 
X=X  COORDINATES  OF  HORZONTAL  GRID  POINTS 
NOTE  -  TAPE1  SHOULD  CONTAIN  STEADY  STATE  INTIAL  CONDITIONS 
OF  BASE  FLOW  FROM  LTRAN2 

TAPE2  SHOULD  CONTAIN  STEADY  STATE  INTIAL  CONDITIONS 
OF  CALIBRATION  FLOW  FROM  LTRAN2 

ON  OUTPUT  TAPE3  WILL  CONTAIN  STEADY  STATE  INTIAL  CONDITIONS 

OF  THE  CURRENT  FLOW  FOR  LTRAN2 

INTIAL  CONDITIONS  FROM  LTRAN2  ARE 

T-TIME 

GLIFT=LIFT 

PU=D I STURBANCE  UELOCITY  POTENTIALS  OF  UPPER  SURFACE 
PL=D I STURBANCE  UELOCITY  POTENTIALS  OF  LOWER  SURFACE 
P=D I STURBANCE  UELOCITY  POTENTIALS  OF  ALL  GRID  POINTS 

READC5. 1 )LMAX» JMAX. JLE. JTE* XSO.  XS1 . EP.  EPO 

1  FORMAT (41 5.  6F10.4) 

|  — XSO 

WRITE(S. 6)LMAX.  JMAX. JLE. JTE. XSO. XS1 » DXS. EP. EPO 

6  FORMAT (/5X.*N0  OF  UER  MESH  POINTS=*. 15. *  NO  OF  HOR  MESH  PTS=*. 15. 

1  *  LE  MESH  PT=*.I5.*TE  MESH  PT=*. 15  .  /5X,*BASE  SHOCK  POSITION^ 

2. F10.4.*CAL  SHOCK  POSITION^.  F10 . 4.  *  SHOCK  DISPLACEMENT^,  FI  0.4/5X. 

3.  *  EP=*»  F10.4* *  EPO=*. F10.4) 

READ(6. 7) CXC I) . 1=1 , JMAX) 

7  FORMAT  C  8E 1 0 . 4 ) 

URITE(6» 11 ) 

11  FORMAT (/5X.*  X  CO-ORDINATES* ) 

URITEC6. 18) CX( I ) . 1=1 . JMAX ) 

18  FORMAT C/5X* 1 OF 12. 6) 

CONS=XSO*d.O-XSO) 

DO  20  1=1. JMAX 
20  TX( I )=0. 0 

DO  30  I=JLE. JTE 

30  TX(I)=(X(I)*(1.0-X(I)))/CONS 
WRITEC6.31) 

31  FORMAT ( /5X» *DISTORSION  C0EF=*) 

WRITE(6.36)(TX(I). 1=1. JMAX) 

36  FORMAT C/5X, 10F12*B) 

REWIND  1 

READ( 1 )T1»LMAX1, JMAX1. GLIFT1.  (PU( J).PL( J). (PCL. J).L=1.LMAX), 

1  J=l» JMAX) 

WRITEC6.2) 

2  FORMAT (/5X.*  DATA  FROM  TAPE1*) 

WRITE (6. 41 )T1. LMAX1. JMAX1. GLIFT1 

41  FORMAT</5X.F10.4.2I5.F10.4) 

DO  10  1=1. JMAX 
WRITEC6.1B)  I* PU( I ).PL( I ) 

16  FORMAT ( /5X . *FOR  COLUMN=*. 15. *PU=*.E15,6* *  PL=*.E15,6) 

WRITE (6. 17) (P( J. I ) » J=l. 10) 

17  FORMAT (/5X. 10E12.4) 

10  CONTINUE 

REWIND  3 

WRITE(3)(PU(I). 1=1. JMAX) 

WRITE(3)CPL(I). 1=1. JMAX) 

DO  40  I=1»LMAX 
40  WRITE(3) CP(I# J). J=l» JMAX) 
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REMIND  2 

READ(2)T2.LMAX2. JMAX2. GLIFTP.  (PUC J) .PLC J). (PCL. J). L=1.LMAX). 
1  J=l. JMAX) 

WRITE(6.  21 ) 

21  FORMAT (/5X»*  DATA  FROM  TAPE  2*) 

WRITE ( G»  4 1 ) T2.  LMAX2*  JMAX2*  GL IFT2 

DO  100  1  =  1*  JMAX 

WRITE (6. 16)I.PU(I)»PL(I) 

100  WRITE<;6.17)(P<J.I>.J=1.10) 

RA=EP/EPO 
REWIND  4 
REWIND  3 
C1=RA*DXS 

READ ( 3 ) ( DUM Y 1  ( I ) ,  I = 1 . JMAX ) 

DO  200  1=1. JMAX 

200  DUMY2(I)=DUMY1(I)*(1.G-C1*TX(I))+RA*(PU(I)-DUMY1(I)* 

1(1. 0~DXS*TX ( I ) ) ) 

WR I TE ( 4 ) (  DUMY2 ( I  )• 1=1* JMAX ) 

READC3) (DUMY1 ( I ) . 1=1 . JMAX) 

DO  300  1=1. JMAX 

300  DUMY2C I )=DUMY1 ( I )*( 1 . 0-Cl*TX( I ) )+RA*(PL( I )-DUMYl (I)* 

1(1. 0-DXS*TXC I ) ) ) 

WR I TE ( 4 ) ( DUM Y2 ( I ) , I = 1 . JMAX ) 

DO  500  1=1. LMAX 

READ ( 3 ) ( DUM Y 1 ( K ) . K= 1 . JMAX ) 

DO  400  J=  .1 .  JMAX 

400  DUMY2( J)=DUMY1( J)*( 1 .0-Cl*TX( J) )+RA*(P( I . J)-DUMYlt J)* 

1(1. 0-DXS*TX( J) ) ) 

500  WRITE(4) (DUMY2(K)*K=1. JMAX) 

T3=T1+RA*(T2-T1) 

GLIFT3=GLIFT1+RA*(GLIFT2-GLIFT1) 

REWIND  4 
REWIND  3 

READC4) (PU( I). 1=1* JMAX) 

READ(4)(PL(I). 1=1. JMAX) 

DO  GOO  1=1. LMAX 
GOO  READ (4) (P( I. J). J=l. JMAX ) 

WRITE ( 3) T3.LMAX1. JMAX 1.GL1FT3.  (PU( J) .PL( J) . (PCL* J) .L=l. LMAX) 
1  J=l. JMAX) 

REWIND  3 

READ(3)T3. LMAX1, JMAX1. GLIFT3*  (PU(J) . PLC J). CP(L. J). L=1.LMAX). 
1  J=l. JMAX) 

WRITE(G.GOl) 

G01  FORMAT (/5X. *  DATA  TO  BE  WRITTEN  ON  TAPE  3*) 

WRITE(6*  41 )T3* LMAX 1 . JMAX 1 . GL I FT3 
DO  700  1=1. JMAX 
WRITE (G. 16) I. PU( I ).PL( I ) 

700  WRITE(G. 17)(P(J.I).J=1*10) 

STOP 

END 
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